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1 Experimental implementation of a two-gene flux match- 
ing system based on negative autoregulation: Materials 
and methods 



1.1 Reactions and domains design 

A graphical sketch of the domain-level design for the self-repression interconnection is shown in 
Figure ST A. The RNA outputs of each genelet are designed so that: 

1) Each RNA output has a domain complementary to its activator strand. 

2) The two RNA species are also complementary. 

These specifications introduce a binding domain between T, and Rj, which is considered another 
off state, as shown in Figure |5T] B. Such a complex is a substrate for RNase H and the RNA 
strand is degraded by the enzyme, releasing the genelet activation domain. We assume that the 
transcription efficiency of an RNA-DNA promoter complex is very low. This hypothesis was not 
experimentally challenged for this specific system; however data shown in Franco et al. j2011|, 
Supplementary Information, show that this assumption is valid for other genelets with the same 
promoter domain. 

The design of a self-inhibiting genelet was first characterized in Kim [2007] . The circuit design 
proposed here, with two-domain RNA transcripts, was originally presented in Franco et al.| (2(308 1 . 

DNA strands were designed by thermodynamic analysis using the Winfree lab DNA design 
toolbox for MATLAB, Nupack [Zadeh et al | [2011] and Mfold |Zuker and Stieglerj |1981| . The 
strands were optimized to yield free energy gains favoring the desired reactions, and to avoid 
unwanted secondary structures and crosstalk. Further constraints on the length and structure 
of the strands, which can affect the transcription efficiency and fidelity, were taken into account 
referring to |Kim| [2007] , Chapter 3.4. 



1.2 Oligonucleotide sequences 

Due to technical constraints of the supplier IDT DNA, Ti — nt and T 2 — nt were shortened with 
respect to the nominal design to have a length of 125 bases. The strands used in the experiments 
are those denoted below as "Short". These modifications did not alter the regulatory domains 
of the transcripts R\ and R 2 - Also the full length of the main transcription products was not 
affected, as verified by gel electrophoresis in Figure S~2 B. 

Ti-nt Full (134-mer) 5'-CTA ATG AAC TAC TAC TAC ACA CTA ATA CGA CTC ACT ATA 
GGG AGA AAC AAG AAC GAC ACT AAT GAA CTA CTA CTA CAC ACC AAC CAC AAC TTT 
ACC TTA ACC TTA CTT ACC ACG GCA GCT GAC AAA GTC AGA AA-3' (not synthesized) 
T^nt Short (125-mer) 5'-Tamra-CT AAT GAA CTA CTA CTA CAC ACT AAT ACG ACT CAC 
TAT AGG GAG AAA CAA GAA CGA CAC TAA TGA ACT ACT ACT ACA CAC CAA CCA CAA 
CTT TAC CTT AAC CTT ACT TAC CAC GGC AGC TGA CAA-3' 

Ti-t (107-mer) 5'-TTT CTG ACT TTG TCA GCT GCC GTG GTA AGT AAG GTT AAG GTA 
AAG TTG TGG TTG GTG TGT AGT AGT AGT TCA TTA GTG TCG TTC TTG TTT CTC 
CCT ATA GTG AGT CG-3' 

A : (35-mer) 5'-TAT TAG TGT GTA GTA GTA GTT CAT TAG TGT CGT TC-3' 
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Figure Si: General reaction scheme representing a transcriptional circuit implementation of the 
two-gene negative feedback scheme for flux matching. Complementary domains have the same 
color. Promoters are in dark gray, terminator hairpin sequences in light gray. The RNA output of 
each genelet is designed to be complementary to its corresponding activator strand. The two RNA 
species are also complementary. A. Desired self-inhibition loops. B. Undesired cross-hybridization 
and RNase H mediated degradation of the RNA-template complexes. 
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T 2 -nt Full (126-mer) 5'-GGT TAA GGT AAA GTT GTG GTT GTA ATA CGA CTC ACT ATA 

GGG AGA AAC AAG TAA GTA AGG TTA AGG TAA AGT TGT GGT TGG TGT GTA GTA 

GTA GTT CAT TAG TGT CGT TCC TGA CAA AGT CAG AAA-3' (not synthesized) 

T 2 -nt Short (126-mer) 5'-TexasRed-GG TTA AGG TAA AGT TGT GGT TGT AAT ACG ACT 

CAC TAT AGG GAG AAA CAA GTA AGT AAG GTT AAG GTA AAG TTG TGG TTG GTG 

TGT AGT AGT AGT TCA TTA GTG TCG TTC CTG ACA AAG TCA GAA-3' 

T 2 -t (99-mer) 5'-TTT CTG ACT TTG TCA GGA ACG ACA CTA ATG AAC TAC TAC TAC 

ACA CCA ACC ACA ACT TTA CCT TAA CCT TAC TTA CTT GTT TCT CCC TAT AGT GAG 

TCG-3' 

A 2 (35-mer) 5'-TAT TAC AAC CAC AAC TTT ACC TTA ACC TTA CTT AC-3' 

Ri (95-mer) 5' - GGG AGA AAC AAG AAC GAC ACU AAU GAA CUA CUA CUA CAC ACC 

AAC CAC AAC UUU ACC UUA ACC UUA CUU ACC ACG GCA GCU GAC AAA GUC AGA AA 

-3' 

R 2 (87-mer) 5'-GGG AGA AAC AAG UAA GUA AGG UUA AGG UAA AGU UGU GGU UGG 
UGU GUA GUA GUA GUU CAU UAG UGU CGU UCC UGA CAA AGU CAG AAA -3' 



1.3 DNA oligonucleotides and enzymes 



All the strands were purchased from Integrated DNA Technologies, Coralville, IA |IDT| 7\ — nt is 
labeled with TAMRA at the 5' end, T 2 —nt is labeled with Texas Red at the 5' end, both activators 
Ai and A 2 are labeled with the IOWA black RQ quencher at the 3' end. The transcription buffer 
mix was prepared prior to each experiment run (two to four samples) using the T7 Megashortscript 
kit (#1354), Ambion, Austin, TX which includes the T7 RNA polymerase enzyme mix, the 
transcription buffer, and rNTPs utilized in the experiments. E. coli RNase H was purchased from 
Ambion (#2292). 



1.4 Transcription protocol 

The templates were annealed with 10% (v/v) lOx transcription buffer from 90°C to 37°C for 1 
h 30 min at a concentration 5-10x the target concentration. The DNA activators were added 
to the annealed templates from a higher concentration stock, in a solution with 10% (v/v), lOx 
transcription buffer, 7.5 mM each NTP, 4% (v/v) T7 RNA polymerase, and .44% (v/v) E. coli 
RNase H. Each transcription experiment for fluorescence spectroscopy was prepared for a total 
target volume of 70 /zl. Samples for gel studies were quenched using a denaturing dye (80% 
formamide, 10 mM EDTA, O.Olg XCFF). 



1.5 Data acquisition and processing 

The fluorescence was measured at 37°C every two minutes with a Horiba/Jobin Yvon Fluorolog 3 
system. Excitation and emission maxima for TAMRA were set to 559 nm and 583 nm, respectively, 
according to the IDT recommendation; for Texas Red the maxima for the spectrum were set to 
598-617 nm. Slit widths were set to 2 nM for excitation and 4 nM for emission. The raw 
fluorescence data $(t) were converted to estimated switch activity by normalizing with respect 
to maximum fluorescence $ max (measured before adding activators and enzymes) and to minimum 
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fluorescence $ min (measured after adding activators and before adding enzymes): 

$(t) - $ ■ 

\ " I mv 



[TiAi](t) = [Tf ot ] • 1 - 



For the adaptation experiments, normalization was done by measuring maximum and minimum 
fluorescence levels at the beginning of the experiment, and assuming that the maximum flu- 
orescence level scales linearly with the change in total fluorescently labeled strands, while the 
minimum is not significantly affected by that variation. We used the formula: 

[TiAi](t) = a[T* ot ] • f 1 



where a is a factor that scales the total amount of template as it varies in the experiment. 

Denaturing polyacrylamide gels (8% 19:1 acrylamide:bis and 7 M urea in TBE buffer, 100 mM 
Tris, 90 mM boric acid, 1 mM EDTA) were run at 67°C for 45 min with 10 V/cm in TBE buffer. 
Samples were loaded using Xylene Cyanol FF dye. For quantitation, denaturing gels were stained 
with SYBR Gold (Molecular Probes, Eugene, OR; #S-11494). In the control lane a 10-base 
DNA ladder (Invitrogen, Carlsbad, CA; #1082-015) was utilized. The DNA ladder 100 bp band 
was used as a control to roughly estimate the concentrations of the RNA species in solution in 



Figure S4 E and F. Gels were scanned using the Molecular Imager FX (Biorad, Hercules, CA) 



and analyzed using the Quantity One software (Biorad, Hercules, CA). 

1.6 Characterization assays 

This section reports experimental results and numerical fits. All experiments were run in trip- 
licates: mean and error bars (standard deviation) are shown in each figure, together with the 
simulated traces (dashed lines) from our fitted model. The full derivation for the model fitted to 
the data is in Section [221 



1.6.1 Genelets in isolation 

Figure 152] A shows the behavior of the two genelets in isolation: we can verify that each genelet 



self-inhibits after the enzymes are added. (For details on the data normalization procedure, 



refer to Section 1.5 ) The concentration of RNA present in solution can be measured through 
gel electrophoresis, as shown in Figure 32 B: lanes 1 and 2 show that free RNA in solution is 
effectively absent. 

1.6.2 Interconnected genelets 

When the two genelets are present in solution in stoichiometric amount, their RNA outputs bind 
quickly to form a double-stranded complex, and therefore the feedback loops become a secondary 
reaction (by design thermodynamically less favorable than the R\ ■ R2 complex formation). As 
shown in Figure[S2]C, the two genelets only moderately self-repress. The total RNA concentration 



in solution is high, as shown in the denaturing gel in Figure [32] B, lanes 3 and 4 



When the templates [T* ot ] and [T^ 1 ] are in different ratios, the system behavior is shown 



in Figure S3 We can plot the resulting initial active template ratio (which corresponds to the 
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0 50 100 150 

Time (min) 



Figure S2: A. Experimental data showing the isolated active genelet concentrations as a function 
of time: the self-inhibition reaction turns the switches off, and the RNA concentration in solution 
is negligible, as verified in the gel electrophoresis data in panel B, lanes 1 and 2 (samples taken at 
steady-state after 2 h). Dashed lines represent numerical trajectories of equations (|5]), using the 
fitted parameters in Table 32 B. Denaturing gel image: lanes 1 and 2 show that the switches in 
isolation self-inhibit and no significant transcription is measured. Lanes 3 and 4 show the total 
RNA amount in samples from the experiment shown at panel C, taken at steady-state after 2 h. 
When the genelets are in stoichiometric amount, their flow rates are already balanced and there 
is only moderate self-inhibition. 
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total template ratio) versus the steady-state one: we find that the system behaves symmetrically 
and the steady-state ratio is close to one across all the initial ratios. Therefore, given open loop 
transcription rates that differ across a factor of 1-3, these results suggest that the system robustly 
matches the flux of R± and i?2- 



1.6.3 Flux adaptation 

If the concentration of [T* ot ] and [A* ot ] is changed over time, the steady-state concentration of 



active genelets adjusts as shown in Figure S4 A and B. Samples from this set of experiments were 



analyzed using a denaturing gel: the results are shown in Figure S4 C and D (corresponding to 



the traces in Figure S4 A and B, respectively) and show the total RNA amount in solution and 
that w as desired (Figure S4 E and F). The RNA concentrations were estimated 



using the DNA ladder as a reference. We are aware that this method may result in inaccurate 
absolute concentration estimates for RNA: however, our objective here is to compare the evolution 
over time of the relative RNA concentrations. Thus, inaccuracy in the determination of the 
absolute amount of RNA produced does not affect the measured outcome of our experiments. 
The adjustment of genelet activity becomes progressively slower over time: the third round of 
adaptation is consistently slower than the previous two. We attribute this slower adaptation to 
various phenomena: 1) Decrease of activity of enzymes over time; 2) Accumulation of incomplete 
degradation products from RNase H hydrolization of RNA in RNA-DNA hybrids: these products 
can be up to 7-8 bases long, and may interfere with the desired inhibition pathways; 3) Abortive 
transcription of RNA, which could also potentially bind to regulatory domains of DNA activators. 
Our hypothesis of accumulation of short products over time is validated by the gels shown in 



Figure S4 C and D (below 60 bases, part of the gel that is not shown, a similar smear is visible). 



1.6.4 Data fitting 

We derived a system of ordinary differential equations (ODEs) starting from mass action kinetics, 



as described in Section 2.2 The ODE system was numerically fitted using MATLAB (The Math- 



Works) to fluorescence data in Figures S2 and S3 Only a subset of the parameters was fit using 



the MATLAB fmincon routine. We fit the total RNA polymerase and RNase H concentrations 
and the rates k TiAi , k TtAiRi} k AiRi , k RlR2 , k RiTj , and the parameters k cat o Nll and k catHtj - This 
specific subset of parameters was chosen because experimental outcomes are chiefly affected by 
branch migration rates (which are tunable by design of the toehold lengths), enzyme concentra- 
tion, and enzyme catalytic rates. The concentration and composition of the transcription enzyme 
mix for the T7 Megashortscript kit are not disclosed by Ambion, but available literature suggests 
that additional enzymes, such as pyrophosphatase, are present in the mix, |Milburn et al.| |U. S.| 
Patent 5256555, 1993 . We neglected reactions associated with the possibly unknown amount 



of pyrophosphatase in the mix. The concentration of RNase H is also not disclosed by Ambion; 
we did not run separate experiments to fit exclusively the degradation rate parameters. A table 
reporting all the parameters is in Section [2] Tables Si and [S2| 
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Figure S3: Concentration of active genelets over time at different total templates concentration. 
The concentration of activators is always stoichiometric to the amount of corresponding template. 
Dashed lines in all the figures correspond to numerical simulations for model ([5]), using the 
parameters in Table |S2 
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Figure S4: A and B. Fluorescent traces showing the adaptation of the active fraction of genelets, 
when the total amount of templates is varied over time. C and D. Samples from the experiments 
shown in panels A and B, respectively, were analyzed with gel electrophoresis. E and F show the 
concentrations of RNA species, estimated from the gel samples. 
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2 Modeling and numerical analysis: two-gene flux match- 
ing system 

2.1 Simple model system: derivation of nullclines and rate matching 
conditions 



We consider a system composed of two generating species Ti 
and T 2 , whose products Ri and R2 interact to form a complex 
P = Ri ■ R 2 . We introduce negative autoregulation to minimize 
the concentration of product that is not used to form the output 
complex (Figure S5). Free molecules of R it i = 1,2, bind to active 
Tj, thereby inactivating it: 



f T1 " 



ION- AY 1 



R1 



Ri 



■tAt* 

T, 



Figure S5: Our two-gene 
negative feedback architec- 
ture 



where T* is an inactive complex. We assume that Tf ot = Ti + T*, and that T* naturally reverts 
to its active state with a first-order rate on. The total amount of Ri is [R\ ot \ = [Ri] + [T*] + [P]. 
The corresponding differential equations are: 

d[T t ] 



dt 
d[Ri 



a i ([Ti] ~ &]) — $i [Ri][Ti], 

Pi [Ti] - k [Ri][Rj] - 5i [Ri][Ti\. 



(1) 



For illustrative purposes, these differential equations are solved numerically. The parameters 
chosen are: <*i = a 2 = 3 • 10" 4 /s, fa = /3 2 = 0.01 /s, 6t = 5 2 = 5 • 10 2 /M/s, and 
k = 2 • 10 3 /M/s. An imbalance in the production rates of Ri and R 2 is created by setting 
[2^(0) = [T{ ot ] = 100 nM and [T 2 ](0) = [T 2 to *] = 200 nM, while [^^(O) = [R 2 ](0) = 0. The 
overall result of this feedback interconnection is that the mismatch in the flow rate of Ri and 
R 2 is reduced, as shown in Figure S6 The flow rate is defined as the derivative of [Rj ot ]- The 
flow rate mismatch is defined as the absolute value of the difference between the two flows. The 
effect of changing the feedback strength, for simplicity chosen as 5i = S 2 , is shown in Figure [57 



the figure shows the mean active fraction of [Xj] and the mean flow mismatch over a trajectory 
simulated for 10 hours. We show the mean over the last hour of the simulation, rather than 
steady-state values. 

It is possible to examine the nullclines relating T\ and T 2 , and find the equilibria Ti and T 2 
as intersection of these nullclines: 

ai(Ti ot -Ti) 



T = 0 



Ri = 0 



Ri 



Ri 



kRj + SiTi 

To simplify the derivation, we set Si = 5 2 = 5, f3± = (3 2 = (3, a± = a 2 = 
expressions for R it we get the following equations (for i = 1,2 and j = 1, 



a. 

2): 



Equating the two 



+ a(Tt ot - - f3T = 0. 
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Figure S6: Numerical simulation showing the solution to the two-gene negative feedback archi- 
tecture for flux matching modeled with equations ([I]). The flow mismatch between Ri and R2 
is shown in the bottom-right panel. 
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Figure S7: Numerical simulation showing genelet activity and mismatch over a range of values 
for the negative feedback parameter 5. 
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/ rptot rp \ 

We can find an expression of the nullclines by introducing a change of variables u = ( 1 T 1 J 

and v = f^zZ^Y and defining fa = fa = (f) 2 k, fa = <xT{ ot , fa = olT\°\ fa = /3T 2 tot , and 
finally ip 3 = $T{ ot : 

u\fav)+u(fav + fa-fa^-)-fa^-=0, (2) 

v 2 (fau) + v(fau + fa-fa-^—)-4> 3 - l — = 0. (3) 

1 + u 1 + u 



The roots of the equations above represent the nullclines of the system. Because all the 
parameters in these equations are positive, there is always a single root. The nullclines are 
numerically solved, for varying 5, in Figure |5"8 



A condition for flow matching at steady-state can be derived as follows: 



R\ — R 2 



0. 



PiT-l — 5iTxR x — f3 2 T 2 — 8 2 T 2 R2- 



Substituting the expressions for Ri and R 2 that can be derived by setting Ti = 0 = T 2 , we get: 

fafx - ai(T*°* - 7\) = (3 2 f 2 - a 2 (T* ot - f 2 ). 
Taking a\ = a 2 = a, (5\ = f3 2 = (5 we get: 



T 2 = T X + (T 2 '°* - T\ ot ) . (4) 

a + p 



The flow matching condition above is shown in Figure S8 orange line (also shown in the main 
paper). If f3 a, i.e., the production of Ri is much faster than the generating species Tj 
inactivation rate, then the condition can be rewritten as: 

f x w f 2 . 



2.2 Differential equations modeling the experimental implementation 

Based on our design specifications and the resulting molecular interactions, we built a model for 
the system starting from the list of occurring chemical reactions. The switches Tj and Tj can have 
three possible states: the on state where activator and template are bound and form the complex 
TiAi, the off state given by free Tf, the off state represented by Rj bound to Tj forming TiRj. 
An off state still allows for RNAP weak binding and transcription. Throughout this derivation, 
the dissociation constants are omitted when assumed to be negligible. It is hypothesized that 
the concentration of enzymes is considerably lower than that of the DNA molecules, allowing the 
classical steady-state assumption for Michaelis-Menten kinetics. 
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Figure S8: Nullclines computed for different values of negative feedback rate 5, and flux matching 
condition (orange) 



Branch migration and hybridization reactions among nucleic acids are, for i E {1,2}, j E {2, 1}: 

Activation 
Inhibition 
Annihilation 
Output formation 
Undesired hybridization 



T 4- A ■ kr i^i T . A- 

Ri + Ti- Ai Ri-Ai + Ti 

Ri "f" Ai A % Ri • Ai 



Ri + Rj A Ri ■ Rj 



h'j ■ V'/'".' Rj • V,. 
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The enzymatic reactions are, for i e {1,2}, j G {2,1}: 

ONii h. 

Transcription: on state RNAP + Tj • A { t- RNAP ■ T t ■ Ai Kcat R Nii RNAP + T { Ai + Rj 

k ONii 

k + 

Transcription: off state RNAP + Tj °J F " RNAP ■ T t kcatc L FF " RNAP + T t + Ri 



^OFFji , 

Transcription: off state RNAP + Rj ■ T t ^ RNAP ■ Rj ■ T { catc L FF ^ RNAP + Rj ■ T { + Rj 

k OFFji 
k+ 

Degradation RNaseH + R, t ■ A, t ^ RNaseH ■ R { ■ Ai - " RNaseH + Ai 

k Hii 
k + .. 

RNaseH + R j ■ T t if RNaseH ■ R j ■ T t Kat "' J1 RNaseH + Tj. 



k 



Hji 



Using the law of mass action, we derive the following ODEs: 

^[T] = - k TiAi [Tj] [A^] + k TiAiRt [Ri] [Ti ■ Ai] - k RjT . [Rj] [Tj] + k catHji [RNaseH ■ R 3 • Tj], 
j\[Ai] = - k TiAi [Ti] [A,] - k AiRi [Ri] [Ai] + k catHii [RNaseH ■ R, ■ A t ], 
d [Ri] = - k RiRj [Ri] [Rj] - k TiAiRt [Ri] [Ti ■ Ai] - k RiTj [Ri] [Tj] - k AiRi [R t ] [A t ] 



dt 

d 



+ k cat oNii [RNAP ■ Ti ■ A,] + KatoFFii [RNAP ■ Ti] + k cat oF Fj i [RNAP ■ Rj • Tj], 
• -Rj] = + k RxRj [Ri] [Rj], 



dt 

^[Rj ■ Tj] = + k RjTi [Rj] [Tj] - k catHji [RNaseH ■ Rj • Tj]. 



(5) 



The molecular complexes appearing at the right-hand side of these equations can be expressed 
using mass conservation: 

[Tj • Ai] = [Tt 1 ] - [Ti] - [Rj ■ Ti], [Ri ■ A t ] = [A**\ - [A t ] - [Tj • A t ]. 

We assume that binding of enzymes to their substrate is faster than the subsequent catalytic 
step, and that the substrate concentration is larger than the total amount of enzyme. These 
assumptions allow us to use the standard Michaelis-Menten quasi-steady-state expressions. The 
Michaelis-Menten coefficients can be immediately defined; for instance, for the ON state of the 
template, define: k MO Nu = k o^+ k ^om l Jhen we fjnd . 

k ONii 



/ + [Ti • A,] + [Tx] + [T 2 • A 2 ] + [T 2 ] + [R 2 -T] + [R, • T 2 ] \ 

kMONll kMOFFU kuON22 KmOFF22 kuOFF2\ kuOFF\2 ) 

[Rx-A x ] t [R 2 -A 2 ] , [R 2 -Ti] , [R!-T 2 y 

kMH\2 



[RNAP tot ] =[RNAP] 
[RNaseH'"} ^RNaseH] (l + ^ + + ^ + 

V KM Hll kMH22 KMH21 



Supplementary Information Appendix - Experiments, Data Processing and Modeling 



16 



We can easily rewrite these equations as [RNAP] = ^ RNA ^ — 1 and [RNaseH] = ^ RNas ^ H 
with a straightforward definition of the coefficients P and H. Finally: 



[RNAP ■ Tt ■ 


A] 


[RNAP tot ] [Ti ■ Ai\ 

" ■ KMONii 


[RNAP ■ Rj 


■Ti\ 


[RNAP tot ] [Rj ■ Ti] 

r ■ KMOFFji 


[RNAP 


■n 


[RNAP tot ] [Ti] 

P ■ k MO FFii 


[RNaseH ■ R t ■ 


M 


[RNaseH tot ] [R { ■ Afi 
H ■ k MHii 


[RNaseH ■ Rj 


■n 


[RNaseH tot ] [Rj ■ 
H ■ kuHji 



which can be substituted in equations (|5]). We note that our numerical fits result in an estimated 
RNAP concentration of about 100 nM: thus, in a subset of our experiments the substrate and 
enzyme concentrations are actually comparable, breaking down one of the assumptions required 
for a quasi-steady-state approximation. Nevertheless, our model overall captures the system 
dynamics satisfactorily. 

The nonlinear set of equations ^ was solved numerically using MATLAB ode23 routine. 

Preliminary numerical analysis Prior to designing DNA strands and testing the system with 
wet lab experiments, we ran numerical simulations using equ ations (|5|) u sing parameters reported 
in Table Si These parameters are consistent with those in Kim et al. 2006] , which were fitted 
from data obtained on a transcriptional system with identical promoter/branch migration design 
specifications and sequence content; thus, we refer the reader to Kim et al. (2006 1 for an accurate 
discussion and comparison to other branch migration, transcription, and degradation parameters 
found in the literature. Figure 59 shows the system trajectories that correspond to zero initial 
conditions for [A] and [R t ], while the complexes [T^] = [T[ ot ] = 100 nM, [T 2 A 2 ] = [T 2 tot ] = 50 
nM, [A\ ot ] = 100 nM and [A^ 1 ] = 50 nM. (The simulation first allows for equilibration of all 
the DNA strands in the absence of enzymes. Only the portion of trajectories after addition of 
enzymes is shown.) The total concentration of enzymes was assumed to be [RN AP tot ] = 80 
nM and [RNaseH tot ] = 8.8 nM, consistent with typical volumes used in our experiments and 



Franco et al 



with enzyme stock concentrations of about 1-1.25 \i M Kim and Winfree 12011 
[201 1| . An example of our numerical simulation results is shown in Figure S9 The behavior of 



the system proved to be consistent with the traces obtained for the simple model system shown 
at Figure [56 



Data fitting results As already indicated in Section 1.6.4 equations ([5]) were fitted to all 
fluorescence data in Figures [S2] and [S3] simultaneously, using MATLAB routine fmincon. Only 



a subset of the parameters was fit: the total RNA polymerase and RNase H concentrations, and 
the rates k TiAi} k TiAl n l} k Al R 1} k RlR , 2 , k RiT] , and the parameters k catONll and k catHlj - Table [52 
lists the results of the data fit; Table [53] reports the constraints used in the fitting procedure. 
Our fits indicate that the hybridization and branch migration rates fitting these experiments are 
higher than what found in Kim et al. |2006| , Franco et al. [201 1| . In particular, the binding rate of 
the RNA species is higher than expected; hybridization rates for complementary RNA strands of 
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Figure S9: Numerical simulation for equations ([5]). Parameters are chosen as in Table Si 

[Tf ot ] = [A\ ot ] = 100 nM, [T* ot ] = [Af*] = 50 nM, [RNAP tot ] = 80 nM, and [RNaseH tot ] = 8.8 
nM. These results are consistent with those of the simple model proposed in equations ([T]), and 



analyzed numerically in Figure S6 
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Table SI: Preliminary Simulation Parameters for Equations ([5]) 



Units: [l/M/s] 


Units: [1/s] 


Units: [M] 


/c Tl ^=4-10 4 


kcatONu = 0.06 


Amcw* = 250 • IO" 9 


k Ti A iRi = 5 • 10 4 


kcatOFFu = 1 • IO" 3 


k MO FF t = 1 • IO" 6 


k AiRi = 5 • 10 4 


kcatOFF,, = -5 • 10~ 3 


k M OF Fl] = 1 • 10~ 6 


k R%T] = 1 • io 4 


kcatHu = -1 


k MHu = 50 • IO" 9 


kjuiii = 1 • 10 6 


kcatHji = -1 


k MHji = 50 • IO" 9 


Units: [M] 


Units: [M] 




[RNAP tot ] = 80 nM 


[i?ATasei? to *] = 8.8 nM 





similar length have (to our knowledge) not been assessed before. The expected concentrations of 
RNA polymerase and RNase H and their k cat values are also higher than in previous studies [Kirn 



et al. |2006| , Franco et al. j2011|, where lower hybridization rates were attributed to the presence 
of incomplete degradation products from RNase H hydrolization of DNA/RNA hybrids. These 
short products, known to have length up to 7-8 bases, may interfere with desired regulatory 
pathways Kim and Winfree |2011| . Because the activity and efficiency of off-the-shelf enzymes 
is known to considerably vary from batch to batch Kim and Winfree [2011| , it is reasonable to 
hypothesize that the RNA polymerase and RNAse H batches used in this set of experiments had 
particularly high activity and low occurrence of incomplete transcription/degradation which can 
slow down other reactions. Indeed, the accumulation of these incomplete products over time may 
be the reason for slower dynamics observed in our adaptation experiments in Figure |5"4 



Table S2: Fitted Parameters for ([5]); other parameters were left unvaried with respect to Table 



SI 



Units: [l/M/s] 


Units: [1/s] 


k T%Ai = 6.6 • 10 5 


kcatONu =0.1, k cat ON 2 2 = 0-09 


k Tl A lRl = 0.7 • 10 5 , k T2A2R2 = 0.6 • 10 5 


kcatHu = - 09 


k Al R! = k A2 R 2 = 4.4 • 10 5 


k C atH 2 i = -03, k cat H 12 = .02 


k RiRj = 4.9 • 10 6 




Units: [M] 


Units: [M] 


[RNAP tot ] = 100 nM 


[RNaseH tot ] = 20 nM 



Supplementary Information Appendix - Experiments, Data Processing and Modeling 



Table S3: Fitting constraints for parameters in Table 



Parameter 


Lower Bound 


Upper Bound 


k T t A, 


10 3 


5 • 10 5 




10 3 


5 • 10 5 




10 3 


5 • 10 5 


^RiRj 


10 3 


5 • 10 5 


kcatONn 


0.01 


0.1 


kcatHn 


0.001 


0.1 


[RNAP tot \ 


15 • 10~ 9 


100 • 10~ 9 


[RNaseH tot ] 


5 • 10" 9 


20 ■ 10" 9 
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3 Modeling and numerical analysis: Scalability of the neg- 
ative feedback scheme for flux regulation 

3.1 Simple model system 

We consider now n generating species T ?; , outputting interacting products R it and we explore dif- 
ferent feedback interconnection topologies. Initial studies on scalability were outlined in |Giordano 



et al. |2013| . ODEs were derived using mass action kinetics and used for numerical simulation 
of three- and four-component networks. Negative autoregulation is implemented, as for smaller 
networks, with a self-repression scheme: when an output is in excess relative to the effectively 
used amount, it down-regulates its own production rate. 

where T* is an inactive complex. We assume that [T- ot ] = [T*] + [T*] and that T* spontaneously 
reverts to its active state with a first-order rate «j. The corresponding differential equation 
describing the template dynamics is the same regardless of the topology: 

Depending on the chosen interaction/binding topology for the products Ri, we find that the 
system exhibits different behaviors, as shown in the following sections. 

3.1.1 Single product topology 

A single product topology occurs when a single complex P is formed by the simultaneous inter- 
action of all the n outputs: 



i=i 

The corresponding differential equations are 



^ = Pi[Ti]-5 i [R i ][T i ]-kf[[R i ], 



dt 

d\P] 
dt 



k niRi 



i=l 



and the total amount of i?j is [Rj ot ] = [Ri] + [T*] + [P]. Figure S10 shows the numerical solutions 
to the ODEs for n = 3 and n — 4. Even though the initial total amounts of Tj are different, 
the concentration of active Xj (bottom left panel) gradually decreases and the flow mismatches 
(namely the differences in absolute value between any two production rates, shown in the bottom 
right panel) are considerably reduced with a fast time response. We can notice that the response 
is slower in the case of 4 interconnected species. The quantity of produced Ri (upper left panel) 
is of course increasing. With respect to the other topologies, as we will see, the single product 
topology leads to a much higher amount of free Ri (upper right panel), which can be considered 
waste because it is not used in the product formation. 
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0 50 100 150 0 50 100 150 



Time (min) Time (min) 

(a) single product, n = 3 




0 50 100 150 0 50 100 150 

Time (min) Time (min) 




0 50 100 150 0 50 100 150 

Time (min) Time (min) 

(b) single product, n = 4 

Figure S10: Example traces from numerical simulations: single product topology, negative feed- 
back scheme. 
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3.1.2 Handshake and neighbor topologies 

A network of n generating species Tj may be designed to produce different subcomponents, 
that may later assemble into a larger product. In this scenario, we can take two extreme cases: 
the neighbor topology, when each output participates in at most two subcomponents, and the 
handshake topology, when each output participates in n — 1 subcomponents. We thus have the 
generation of pairwise products P^; in the handshake case i,j = l,...,n, j ^ i, while in the 
neighbor case i — 1, n, j — i — l,i + 1 and when i — 1, % — 1 = n, when i — n, i + 1 = 1, 
to close the loop. It is worth noticing that, in the case n = 3, the two topologies coincide. The 
reactions corresponding to product generation are 

Ri ~\~ Rj Pij , 

which lead to the following ODEs: 

= a [Tj] - 5 t mm - J2 hj mm, 



dt 



d[Pij 



dt 



shows the numerical 



The total amount of Ri is [Rf 1 ] = [R,] + [T*\ + J2j[ p ij}- Fi g ure Sl1 
solutions to the ODEs for n = 3, and for n = 4 in the handshake connection case. As for the 
single product topology, even though we initially have different total amounts of active T it the 
concentration of active Tj decreases and the flux mismatches are considerably reduced with a fast 
time response. Although the quantity of produced R, L is increasing, the feedback control reduces 
and keeps bounded the amount of free R it which can be considered waste. 

3.1.3 Parameters 

The parameters chosen in our simulations are: = 2 • 10 3 /M/s for the handshake/neighbor 
topology and k = 6-10 3 /M/s for the single product topology, Si = 5- 10 3 /M/s, «j = 3 - 1CT 4 /s, 
Pi = 1 ■ 1(T 2 /s, [Tl ot ] = 100 nM, [T 2 * ot ] = 200 nM, [T*°'] = 300 nM, [T{ ot ] = 150 nM. An 
imbalance in the production rates of Ri is created by setting [Tj](0) = [T/°*], while [-R;](0) = 0. 

3.1.4 Performance overview of the different topologies as a function of key parame- 
ters 

We numerically explored the behavior of the different network topologies for n = 4 as a function 
of the feedback parameter 5 and of the rate of activation a. Figures [Sl2| |513| and |514| show the 



network response in terms of active percentage of Tj ([Tj]/[Tj iot ] ■ 100), flow mismatch (computed 
as in the previous cases) and response time (defined as the time it takes for the active T trajectory 
to go from [Tj(0)] - 10%A to [T(0)] - 90%A, where A is the difference between its initial value 
[T(0)] and its steady state value). We solved the differential equations for a time span of 10 
hours and averaged the trajectories for active T and for the computed mismatch over the last 
simulation hour. 5 varies logarithmically from a tenth to a thousand times its nominal value; a 
varies from a hundredth to five times its nominal value. In each figure, pink squares mark the 



nominal behavior of the system (all parameters are identical to those listed in Section 3.1.3) 
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0 50 100 150 0 50 100 150 



Time (min) Time (min) 

(a) handshake/neighbor, n = 3 




0 50 100 150 0 50 100 150 

Time (min) Time (min) 




0 50 100 150 0 50 100 150 

Time (min) Time (min) 

(b) handshake, n = 4 

FigureSll: Example traces from numerical simulations: handshake/neighbor topologies, negative 
feedback scheme. 
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In all network topologies, a large negative feedback parameter 5 yields a lower mismatch and 
decreases the response time; however, large 5 clearly reduces the steady state activity of Tj. In 
the handshake and neighbor topologies, a larger value of the spontaneous reactivation parameter 
a yields higher T steady state activity, a larger mismatch, and a shorter response time. On the 
contrary, in the single product topology larger a, despite yielding higher T steady state activity, 
dramatically increases the response time, while the mismatch does not monotonically increase. 
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Figure S14: Simulations for the negative feedback, neighbor topology: parameter sensitivity analysis. 
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4 Positive feedback architecture for a two-gene system. 
Modeling and a viable experimental implementation 



4.1 



Simple model system: derivation of nullclines and rate matching 
conditions 



As done for the negative feedback architecture, we consider a 
system composed of two generating species T\ and T 2 , whose 
products Ri and R 2 interact to form a complex P = R\ ■ R 2 . We 
devise a positive feedback interconnection where product in excess 



upregulates the product in shortage (Figure S15). Free (and thus, 
in excess) molecules of Ri bind to inactive Tj and activate it: 

Ri + Tj Tj 
2* 2 77, 




Figure S15: Our two-gene 
positive feedback architec- 
ture 



where again T* is an inactive complex and [T/ ot ] = [Tj] + [T*]. The total amount of Ri is 
[Rj ot ] = [Ri] + [Tj] + [P]. We now assume that Tj naturally reverts to its inactive state with rate 
«j. The corresponding differential equations are 



d[Tj] 

dt 
d[Rj] 
dt 



[T]+^ [iycm-ro, 
hin-kmRA-Sijim^-im. 



(6) 



This system was initially considered in Franco |2012 . The above differential equations were 
solved numerically. The parameters were chosen for illustrative purposes as oi\ = a 2 = 3-10 -4 /s, 
fa = f3 2 = 0.01 /s, 5 t = 5 2 = 5 • 10 2 /M /s, and k = 2- 10 3 /M /s. The total amount of templates 
was chosen as [Tf 5 *] = 100 nM, [T 2 * ot ] = 200 nM. The initial conditions of active [Tj] are set as 
[Ti](0) = 10 nM and [T 2 ](0) = 160 nM, while [Ri](0) = [i2 2 ](0) = 0. Example traces are shown 
in Figure S16| (a modified version of this figure is also in the main paper). Each product's flux 
rate is defined again as the derivative of [R\ ot ]- The flux mismatch is defined as the absolute 
value of the difference between the two flux rates. The effect of changing the feedback strength, 
where for simplicity 5\ = 5 2 , is shown in Figure [Sl6] B and C, which plots the active fraction 



of [Tj] and the flux mismatch averaged over the last hour of a 10 hours simulation. The right 
panel in Figure 516 seems to indicate that the flux mismatch of the two circuits is minimized for 
a certain range of 5 around the nominal value of 5 = 5 • 10 2 . 

The nullclines of the system in the Ti~T 2 space can be calculated as done for the negative 
feedback design. Taking equations (|6]), we find: 
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10 10 

6 [/M/s] 



10 10 

6 [/M/s] 



Figure S16: A: Example numerical simulation showing the time evolution of the source species in 
the positive feedback architecture (Figure S15 modeled with equations ([6]). B: Active concentra- 
tion of source species as a function of the positive feedback rate 5. C: Flow mismatch between 
Ri and R2 as a function of 5. Dark circles indicate the value of 5 used in panel A. 
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To simplify the derivation, we set <5 12 = 6 2 \ = 5, fi\ = p~2 — ft, «i = «2 = ex. Equating the two 
expressions for R it we get the following equations (for i — 1, 2 and j = 1,2): 



a \ 2 k ( Ti \ (_JE> 
5 J \Tf ot -Tj {t^-Tj 



* UsrV UisrSr +^-^ = 0- ( 7 ) 



We can find an expression of the nullclines by introducing a change of variables z = ^ Tto Ti Ti j 
and w = (533^). and defining fa = fa = (f) 2 *;, 0 2 = alf 4 , ^2 = «T 2 tot , 0 3 = /3T 2 toi , and 



finally ^ 3 = /3T* ot 

z 2 f<^)i?; N ) 4- z(StW 4- (Ao — fa 

1 + w l + w 



Z 2 {fav) + Z {faw + 0 2 - fa — )-03— = 0, (8) 



w\faz)+w(faz + fa-^-^—)-^-^— =0. (9) 



The roots of the equations above represent the nullclines of the system. Because all the 
parameters in these equations are positive, there is always a single root. The nullclines are 
numerically solved, for varying 5, in Figure S17 A condition for flow matching at steady-state 
can be derived as follows: 



R\ — R'2 = 0, 

0 1 T 1 - faR^T™ - T 2 ) = (3 2 T 2 - 5 12 R 2 {T{ ot - T x 



Substituting the expressions for Ri and R 2 that can be derived by setting T\ = 0 = f 2 , we get: 

$21 - - $12 

PiTx - ^a 2 T 2 = /3 2 T 2 - —a^. 

0l2 $21 

Taking a x = a 2 = a, f} x = p 2 = /3, and 5 12 — 5 2 i = 6 we get: 



f 2 = f 1 . 



(10) 



This flow matching condition is shown in Figure 517 in the red dashed line. Decreasing a 
(inactivation rate for the generating species) or increasing 5 (speed of the positive feedback), 
with respect to the nominal values chosen here, causes the equilibrium of the system to be 
pushed toward the upper right corner of Figure 517 Moreover, when decreasing a or increasing 
5 the system reaches equilibrium on a timescale in the order of several dozens of hours. Explicit 
tradeoffs on the effects of a and 5 may be found by further analysis on the nullclines and on the 
locus of equilibria in equation ([7]). 



4.2 A possible experimental implementation of a two-gene positive 
feedback scheme 

The experimental implementation of our positive feedback scheme using transcriptional networks 
presents several challenges. Here we present its general idea. A viable strand design scheme is in 
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Nullclines for varying 8 



200 



150 



CM 



100 




60 80 
T1 [nM] 

Figure S17: Numerical simulation: nullclines of the positive feedback scheme (|6]) in the 7\- 
T 2 plane, calculated for different values of S finding the roots of equations (|8|) and (|9~]). The 



equilibrium corresponding to the set of nominal parameters (trajectories in Figure S16 A) is 



circled in black. The flow matching condition (10) is shown in the orange line. The flow matching 
condition is satisfied by the equilibria T\ and T 2 for 5 = 5 • 10 3 . 
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Figure S18 A. Both genelets are constitutively inhibited by a DNA inhibitor Jj. Each RNA output 
Ri is designed to bind to the inhibitor Ij (domains indicated as qj-cij-tj), thereby releasing the 
activator Aj for binding to Tj. Because Ri should also cover the active domain of Rj in the 
formation of P, then Ri must also be complementary to Ai (domains t^-a^-q'A: therefore, this 
design is structurally affected by binding of RNA to templates (as for the self-repressing circuit), 



and by RNA-mediated self-inhibition loops, as shown in the reaction scheme in Figure 518| C 



The entity of these design pitfalls depends on the length and sequences of the complementarity 
domains shared by Ri and Rj. For instance, we could avoid inserting in the RNA species the 
toehold sequences t\, t\, t 2 , and t' 2 to minimize the self inhibition; however, this would facilitate 
the formation of complexes A4 • Jj • Rj that would slow down the release of A^. 

Preliminary experiments on this design, reported in Franco |2012 , show that the issues de- 



scribed above are significant. In particular, the design could be improved if the self-inhibition 
pathways were minimized: this was attempted, without conclusive success, by increasing the 
concentration of DNA inhibitors, the concentration of RNase H, and by lengthening the length 
of toeholds for Ai and U. Experiments also highlighted the possibility of "leaky" transcription of 
inhibited switches. We refer the reader to Franco [2012 , Chapter 1, for further details. Here, 



we only describe our numerical analysis, which suggests that the scheme has the ability to match 
transcription rates of two cross-activating genelets when we choose plausible reaction parameters. 

4.2.1 Modeling 

To construct a dynamic model for the cross-activating circuit represented in Figure 518 A, we 
start from a list of all the chemical reactions that can occur, 



Activation 


T^A^T^Ai 


Inhibition 


Ti-Ai + h kT ^ h Ti + h- Ai 


Annihilation 


Ai -\- Ii -4 1 Ai ■ Ii 


Release 


Ri + Aj ■ Ij kR ^ I] Ri ■ Ij + Ai 


Annihilation 


Ri • //" / />W, 


Output formation 


Ri + Rj A J Ri ■ Rj 


Undesired interactions 


Ri + Ai A 1 Ri ■ Ai 




Ri + Tj k % T i Ri ■ Tj 
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p q2 a2 t2 t1' a1' q1' _ 

= q 2' a2' t2' t1 a1 q1 



R1 q2 a2 t2 t1' a1' ql* 

1 rnap T1 on 



ql a1 t1 t2' a2' q2' - R2 

a2 . ', rnap T2 on 








t2 a2 



Figure S18: General reaction scheme of the transcriptional circuits implementation for the positive 
feedback scheme in Figure [Sl5| Complementary domains are represented with the same color. 



Promoters are colored in dark gray, while hairpin terminator sequences are in light gray. A. Desired 
cross-activation loops. The activation reaction arrows are colored in red. B. Undesired cross- 
activation and RNase H-mediated degradation of the RNA-template complexes. C. Undesired 
self-inhibition. The inhibition pathway in cyan arrows nominally should not occur, since there 
is no exposed toehold to favor it. However, this reaction has been observed in preliminary 
experiments not shown in this manuscript and is therefore also included in the models. 



Supplementary Information Appendix - Experiments, Data Processing and Modeling 



32 



Transcription: on state RNAP + Ti-Ai Z RNAP ■ T { ■ A t cat Z " RNAP + T t -Ai + R { 

k ON ii 

k+ 

OFF,; , 

Transcription: off state RNAP + % Z RN AP ■ Ti cat ° FF * RN AP + Ti + Ri 

k~ 

OFF^ 

k+ 

OFF i:j , 

RNAP + Ri ■ Tj Z RNAP-Ri-Tj cat ^ Fii RNAP + R { ■ Tj + R, 

k OFF i3 
k+ 

I j u 

Degradation RNaseH + R, t ■ Ij Z RNaseH ■ R { ■ Ij RNaseH + Ij 

3 

k' 



kratH A . 



H A . 

RNaseH + R { ■ A, Z RNaseH ■ Ri ■ A { RNaseH + A { 

k: 



k TT 



RNaseH + Ri ■ Tj Z RNaseH ■ R { ■ Tj ^ ^ RNaseH + Tj. 



The resulting set of ordinary differential equations is: 

-[Ti\ = ~ k T t A, Ft] [M - k R . Ti [Rj] [Ti\ + Hmu [Ti • M [h] + k catHn [RNaseH ■ R 3 • TJ, 
j t [A t ] = - k TiAi [Ti] [Ail ~ W< [Ail [hi ~ k RlAl [Ri] [A^ + k cat H Ai [RNaseH ■ R { ■ A<], 

j t [Ii] = - k AiIi [A,] [hi - k TtAih [Ti ■ Ai] fa] - k RjIi [Rj] [Ii] + k catHh [RNaseH ■ Rj ■ 
j t m = - k RiAjI . [Ri] [A, ■ Ij] - k RiRj [R t ] [Rj] - k RiTj [Ril [Tj] - k RiIj [R t ] [/,] - k RiAi [Ri 

d 



+ k catONm [RNAP ■ Ti ■ A,] + k catOFFi [RNAP ■ Tj + k catOFF]1 [RNAP ■ Rj ■ Tj] 
[Ri ■ T^ = + k RiTj [Ri] [Tj] - k catHT . [RNaseH ■ R 4 ■ Tj], 



dt 

d [Ri-R j l= + k RlRj [Ril [Rj]. 



dt 



(11) 



As previously done for the self-inhibiting circuit model, we can express the enzyme-substrate 
complexes using the Michaelis-Menten approximation. For the RNAP substrate, for instance, we 
find: 

[RNAP • * • ^ = / 1|r dlS i [Ri-Tj] \ ' d2) 

V k MON .. ' k MOFF . ~ k MOFF .. 

Analogous expressions can be derived for all other complexes. 



Equations (11) are numerically solved using the MATLAB ode23s solver. Table S4 shows the 



parameters used for the simulations. Such generic parameters are consistent with those in Kim 
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et al. j2006|. For simplicity we assume that the two sub-circuits are symmetric and have the 
same binding rates. We check the behavior of the system by creating an imbalance in the total 
concentration of inhibitors: [T*°<] = [A\ ot ] = 50 nM, [T* ot ] = [A* ot ] = 100 nM, while [I\ ot } = 20 
nM and [Itf 1 ] = 120 nM. The simulation first allows for equilibration of all the DNA strands in 
the absence of enzymes. The plot shows the trajectories after addition of the enzymes, whose 
total concentration is assumed to be [RNAP tot ] = 80 nM and [RNaseH tot ] = 8.8 nM, based on 
typical experimental conditions. As noted before for the self-inhibitory scheme, the concentration 
of RNAP is not negligible relative to the total amount of genelets present and therefore the 
Michaelis-Menten approximation may not be accurate in this case. The simulation results are 
shown in Figure [Sl9] and are consistent with the traces obtained for the simple model system 
shown at Figure 516 A: the templates cross-activate and reach an equilibrium where the flow of 
total RNA is matched. A comparison between the performance of the transcriptional negative 



and positive feedback circuits models was also done in Franco and Murray [2008 1 



Table S4: Parameters for the Initial Numerical Analysis of the Cross Activating Circuit 



Units: [l/M/s] Units: [1/s] Units: [M] 



kr t A t = 


4- 10 4 


kcatONu = 


= 0.06 


kMON u = 


= 250 • 10" 9 




= 5 • 10 4 


kcatOFFi 


= 1 • 10~ 3 


kyiOFFi 


= 1 • 10~ 6 


kAili = 


5- 10 4 


kcatOFFij 


= 1 • 10" 3 


kyiOFFij 


= 1 • 10" 6 


k RjA,h 


= 5 • 10 5 


kcatHj. = 


0.1 


kuH^ = 


50 • 10~ 9 


k R>h = 


5 • 10 5 


kcatH T . — 


0.1 


kuH T . = 


50 • 10" 9 


k RiTj = 


1 • 10 3 


kcatH A . = 


0.1 


kMH At = 


50 • 10- 9 


kRiA x = 


1 • 10 3 










k R l R ] = 


= 2 • 10 5 











5 Numerical scalability analysis of our simplified positive 
feedback scheme model for flux regulation 

Here we report the mathematical models for positive feedback topologies in the case of n gener- 
ating species Tj. This numerical study was initially outlined in Giordano et al. |2013 . The ODE 
systems were derived using mass action kinetics and used for the numerical simulation of the 
proposed topologies in the case of three-component networks. Positive feedback is implemented, 
as for smaller networks, with a cross-activation scheme: when an output is in excess (not used 
in the product formation), it increases the generation rate of all the other outputs it forms a 
product with: 
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Figure S19: Numerical simulation for equations (11). Parameters are chosen as in Table S4 

\T[ ot ] = [A\ ot ] = 50 nM, [T 2 tot ] = [A%*] = 100 nM, while [I{ ot ] = 20 nM, and [1^] = 120 
nM. [RNAP tot ] = 80 nM and [RNaseH tot ] = 8.8 nM. These numerical results are in general 
consistent with those obtained for the simple model ([6]), shown in Figure S16 A. 



where T* is an inactive complex. We assume that [T- ot ] = [Tj] + [T*] and that the active complex 
Tj naturally inactivates with a first order rate a.^ 

5.1 Single product topology 

In a single product topology, a single complex P is concurrently formed by all the n outputs: 

n 



i=l 



The corresponding differential equations are 

dpi] 



dt 



-on 



^rr- = a n - k fim - s> m HdT^] - [Tj]) 



dt 

d[p] 

dt 



(13) 



i=l 



k l[{R t ] 



and the total amount of Ri is [Rf 1 ] = [Ri] + X^yjPj] + [P]- The simulation results, in Figure 
S20 (a), show that also this feedback strategy is effective. The concentrations of active % 
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asymptotically decrease and flow mismatches reduce, yet the time response is slower than in the 
negative feedback case. With respect to negative feedback, there is also a higher Ri production. 

5.2 Handshake and neighbor topologies 



Subcomponents generation is expressed by the reaction Ri + Rj " and positive feedback acts 
on gene i due to gene j: Ri + T* ^ Tj. The differential equations are 

^§ = -a t [T i ] + y £6 ij [R j ]{pt*\-[Ti]) 
j 

^ = a m - £ h 3 miRj] -i> [^wf] - \m ( l4 ) 



dt 
dt 



kij [Ri] [Rj 



and the total amount of Ri is [Rj ot ] = [Ri]+J2j[ T j\+J2j[Pij]- We remind that the two topologies 



coincide in the case n = 3. The simulation results are shown in Figure S20 (b). The concentration 
of active genes decreases and the flux mismatches are reduced, but the response time is still longer 
than in the negative feedback architecture. Moreover, there is a higher Ri production than in 
the negative feedback case. We can note that the handshake/neighbor connection generates less 
waste (unused Ri) than the single product interconnection. 



5.3 Parameters 

For the numerical solution, the parameters chosen are: = 2 • 10 3 /M/s for the hand- 
shake/neighbor topology and k = 6 • 10 3 /M/s for the single product topology, 5ij = 50 /M/s, 
ai = 3 • 1(T 4 /s, ^ = 1 • 1(T 2 /s, [T{ ot ] = 100 nM, [T* ot ] = 200 nM, [T 3 to *] = 300 nM. An 
imbalance in the production rates of Ri is created by setting [Tj](0) = [T<; ot ], while [i?t](0) = 0. 



5.4 Performance overview of the different topologies as a function of 
key parameters 



Using Figures 522| and S23 as a support, we can compare the performance of the positive feedback 
strategy for networks with n — 3. These topologies are shown in Figure S21 for n = 3 the 
handshake and neighbor topology coincide, Figure [S2T] B. 



We numerically analyzed the network response in terms of active percentage of T i: mean 
flow mismatch and response time, defined as previously done for negative feedback topologies. 
We solved the differential equations for a time span of 10 hours and averaged the trajectories 
for active T, t and for the computed mismatch over the last simulation hour. We examined the 
sensitivity to variations in S, the feedback strength, and in a, the rate of spontaneous inactivation 
of Tf. § varies from a hundredth to a hundred times its nominal value; a varies from a hundredth 
of its nominal value to twice its nominal value, and up to five times its nominal value in the 
response time analysis. In each figure, a pink square highlights the system behavior when the 



nominal parameters in Section 5.3 are adopted 
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(b) handshake/neighbor topologies, n = 3 

Figure S20: Example traces from numerical simulations: positive feedback scheme. 
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Figure S21: A: Single product topology. B: Handshake/neighbor interconnection. 
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Figure S22: Positive feedback, single product topology: parameter sensitivity analysis. 
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In all network topologies, an increase in the spontaneous inactivation parameter a yields a 
lower mismatch, decreases the response time and considerably reduces the steady state activity of 
Tj. In the handshake/neighbor topology, an increase in the positive feedback parameter 5 yields a 
significantly higher Tj steady state activity and a larger mismatch; in the single product topology, 
instead, the steady state activity of T is quite low and almost insensitive to variations in S and 
the mismatch is almost independent of 5. When 5 increases, the response time decreases in the 
single product topology, while it does not have a monotone behavior in the handshake/neighbor 
topology. 
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